home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Meeting Pearls 1
/
Meeting Pearls Vol 1 (1994).iso
/
amok98-106
/
amok101
/
matrix
/
matrix.mod
< prev
next >
Wrap
Text File
|
1994-02-08
|
9KB
|
437 lines
(*--------------------------------------------------------------------
:Program. MATRIX.mod
:Author. Michael Meyer
:Address. Würmseestrasse 28, 81476 München
:EMail. ZNet: M.MEYER@AMC.ZER.SUB.ORG
:EMail. Fido: Michael Meyer 2:2480/21.100
:Version. 1.0
:Date. 6 Feb 1994
:Copyright. 1994 Michael Meyer
:Contents. Dieses Modul unterstüzt das Rechnen mit Matrizen
:Language. Oberon-2
:Translator. Amiga Oberon 3.00d
--------------------------------------------------------------------*)
(* Matrix = elements( [0, 0] [0, 1] [0, 2] ... [0, n-1]
* [1, 0] [1, 1] [1, 2] ... [1, n-1]
* [2, 0] [2, 1] [2, 2] ... [2, n-1]
* . . . .
* . . . .
* . . . .
* [m-1,0] [m-1,1] [m-1,2] ... [m-1,n-1] )
*)
MODULE MATRIX;
IMPORT BT * := BasicTypes;
TYPE
MATRIX * = POINTER TO MATRIXDesc;
MATRIXDesc * = RECORD (BT.GROUPDesc)
m, n: LONGINT;
elements: POINTER TO ARRAY OF ARRAY OF LONGREAL;
END;
PROCEDURE Create * (m, n: LONGINT): MATRIX;
(*
* require
* m>0 & n>0
*)
VAR
new: MATRIX;
i, j: LONGINT;
BEGIN
NEW(new);
NEW(new.elements, m, n);
new.m := m;
new.n := n;
FOR i := 0 TO m-1 DO
FOR j := 0 TO n-1 DO
new.elements[j, i] := 0.0;
END;
END;
RETURN new;
END Create;
PROCEDURE NeutralMatrix * (m: LONGINT): MATRIX;
(*
* require
* m>0 & n>0
* a.DimM() = a.DimN()
*)
VAR
new: MATRIX;
i, j: LONGINT;
BEGIN
NEW(new);
NEW(new.elements, m, m);
new.m := m;
new.n := m;
FOR i := 0 TO m-1 DO
FOR j := 0 TO m-1 DO
new.elements[j, i] := 0.0;
END;
END;
FOR i := 0 TO m-1 DO
new.elements[i, i] := 1.0;
END;
RETURN new;
END NeutralMatrix;
PROCEDURE (a: MATRIX) Get * (m, n: LONGINT): LONGREAL;
(*
* require
* a.DimM()>m &
* a.DimN()>n;
* m>0 &
* n>0
*)
BEGIN
RETURN a.elements[m, n];
END Get;
PROCEDURE (a: MATRIX) Put * (m, n: LONGINT; x: LONGREAL);
(*
* require
* a.DimM()>m &
* a.DimN()>n;
* m>0 &
* n>0
*)
BEGIN
a.elements[m, n] := x;
END Put;
PROCEDURE (a: MATRIX) DimM * (): LONGINT;
BEGIN
RETURN a.m;
END DimM;
PROCEDURE (a: MATRIX) DimN * (): LONGINT;
BEGIN
RETURN a.n;
END DimN;
PROCEDURE (a: MATRIX) Add * (b: BT.GROUP): BT.GROUP;
(*
* require
* a.DimM() = b(MATRIX).DimM()
* a.DimN() = b(MATRIX).DimN()
*)
VAR
new: MATRIX;
i, j: LONGINT;
BEGIN
WITH b: MATRIX DO
new := Create(a.m, a.n);
FOR i := 0 TO a.m-1 DO
FOR j := 0 TO a.n-1 DO
new.elements[i, j] := a.elements[i, j] + b.elements[i, j]
END;
END;
END;
RETURN new;
END Add;
PROCEDURE (a: MATRIX) Sub * (b: BT.GROUP): BT.GROUP;
(*
* require
* a.DimM() = b(MATRIX).DimM()
* a.DimN() = b(MATRIX).DimN()
*)
VAR
new: MATRIX;
i, j: LONGINT;
BEGIN
WITH b: MATRIX DO
new := Create(a.m, a.n);
FOR i := 0 TO a.m-1 DO
FOR j := 0 TO a.n-1 DO
new.elements[i, j] := a.elements[i, j] - b.elements[i, j]
END;
END;
END;
RETURN new;
END Sub;
PROCEDURE (a: MATRIX) Mul * (b: BT.GROUP): BT.GROUP;
(*
* c = a * b
*
* require
* a.DimN() = b(MATRIX).DimM()
*
* => c.DimM() = a.DimM()
* c.DimN() = b(MATRIX).DimN()
*)
VAR
new: MATRIX;
i, j, o: LONGINT;
c: LONGREAL;
BEGIN
WITH b: MATRIX DO
new := Create(a.m, b.n);
FOR i := 0 TO a.m-1 DO
FOR j := 0 TO b.n-1 DO
c := 0.0;
FOR o := 0 TO a.n-1 DO
c := c + a.elements[i, o] * b.elements[o, j];
END;
new.elements[i, j] := c;
END;
END;
END;
RETURN new;
END Mul;
PROCEDURE (a: MATRIX) SkalarMul * (r: LONGREAL): BT.GROUP;
(*
* multiply each element of a by r.
*)
VAR
new: MATRIX;
i, j: LONGINT;
BEGIN
WITH a: MATRIX DO
new := Create(a.m, a.n);
FOR i := 0 TO a.m-1 DO
FOR j := 0 TO a.n-1 DO
new.elements[i, j] := r * a.elements[i, j];
END;
END;
END;
RETURN new;
END SkalarMul;
PROCEDURE (a: MATRIX) isEqual * (b: BT.COMPAREABLE): BOOLEAN;
VAR
i, j: LONGINT;
BEGIN
WITH b: MATRIX DO
IF ~((a.m = b.m) & (a.n = b.n)) THEN RETURN FALSE END;
FOR i := 0 TO a.m-1 DO
FOR j := 0 TO a.n-1 DO
IF a.elements[i, j] # b.elements[i, j] THEN RETURN FALSE END;
END;
END;
END;
RETURN TRUE;
END isEqual;
PROCEDURE (a: MATRIX) isDimEqual * (b: BT.COMPAREABLE): BOOLEAN;
BEGIN
WITH b: MATRIX DO
IF ~((a.m = b.m) & (a.n = b.n)) THEN
RETURN FALSE;
ELSE
RETURN TRUE;
END;
END;
END isDimEqual;
PROCEDURE (a: MATRIX) TransposeMatrix * (): BT.GROUP;
VAR
i, j: LONGINT;
new: MATRIX;
BEGIN
new := Create(a.n, a.m);
FOR i := 0 TO a.m-1 DO
FOR j := 0 TO a.n-1 DO
new.elements[j, i] := a.elements[i, j];
END;
END;
RETURN new;
END TransposeMatrix;
PROCEDURE FindPivot (a: MATRIX; c, l : LONGINT): LONGINT;
(*
* find Pivot-element
* in column c
* from line l on
*
* returns line of Pivot-element
*)
VAR
i : LONGINT;
pivot, next : LONGREAL;
line : LONGINT;
BEGIN
pivot := ABS(a.elements[l, c]);
line := l;
FOR i := l+1 TO a.m-1 DO
next := ABS(a.elements[i, c]);
IF next > pivot THEN
pivot := next;
line := i;
END;
END;
RETURN line;
END FindPivot;
PROCEDURE ChangeTwoLines (a: MATRIX; l1, l2, c: LONGINT);
(*
* change two lines
* from column c to column a.DimN()-1
*)
VAR
i : LONGINT;
x : LONGREAL;
BEGIN
FOR i := c TO a.n-1 DO
x := a.elements[l1, i];
a.elements[l1, i] := a.elements[l2, i];
a.elements[l2, i] := x;
END;
END ChangeTwoLines;
PROCEDURE MulLineAddLine(a: MATRIX; l1, l2, c: LONGINT; s: LONGREAL);
(*
* multiplicate each element of line l1 from column c with
* skalar s and add the result to line l2
*)
VAR
i : LONGINT;
x : LONGREAL;
BEGIN
FOR i := c TO a.n-1 DO
a.elements[l2, i] := a.elements[l2, i] + a.elements[l1, i] * s;
END;
END MulLineAddLine;
PROCEDURE (a: MATRIX) Neg * (): BT.GROUP;
(*
* require
* a.DimM() = a.DimN()
* a.Det() # 0
*)
VAR
new, newN, return: MATRIX;
i, j, k: LONGINT;
ln: LONGINT;
t1, t2, skalar, pivot, det: LONGREAL;
sign: BOOLEAN;
BEGIN
new := a;
newN := NeutralMatrix(a.m);
return := Create(a.m, a.n);
(*------------------ Gauss ------------------*)
sign := TRUE;
FOR i := 0 TO new.m-2 DO
ln := FindPivot(new, i, i);
pivot := -1.0 * new.elements[ln, i];
IF ~(ln = i) THEN
sign := ~sign;
ChangeTwoLines(new, i, ln, i);
ChangeTwoLines(newN, i, ln, 0);
END;
FOR j := i+1 TO new.m-1 DO
skalar := new.elements[j, i]/pivot;
MulLineAddLine(new, i, j, i, skalar);
MulLineAddLine(newN, i, j, 0, skalar);
END;
END;
IF ~sign THEN
FOR i := 0 TO new.n-1 DO
new.elements[new.m-1, i] := -1.0 * new.elements[new.m-1, i];
newN.elements[new.m-1, i] := -1.0 * newN.elements[new.m-1, i];
END;
END;
(*------------------ Jordan ------------------*)
FOR i := new.m-1 TO 1 BY -1 DO
t1 := -1.0 * new.elements[i, i];
FOR j := i-1 TO 0 BY -1 DO
t2 := new.elements[j, i] / t1;
MulLineAddLine(new, i, j, j, t2);
MulLineAddLine(newN, i, j, 0, t2);
END;
END;
FOR i := 0 TO new.m-1 DO
FOR j := 0 TO new.n-1 DO
return.elements[i, j] := newN.elements[i, j]/new.elements[i, i];
END;
END;
RETURN return;
END Neg;
PROCEDURE (a: MATRIX) Gauss * (): BT.GROUP;
(*
* require
* a.DimM() = a.DimN()
*)
VAR
new: MATRIX;
sign: BOOLEAN; (* sign of the matrix: TRUE=plus and FALSE=minus *)
i, j: LONGINT;
ln: LONGINT; (* line-number in which Pivot-element is *)
skalar, pivot: LONGREAL;
BEGIN
new := Create(a.m, a.n);
new := a;
sign := TRUE;
FOR i := 0 TO new.m-2 DO
ln := FindPivot(new, i, i);
pivot := -1.0 * new.elements[ln, i];
IF ~(ln = i) THEN
sign := ~sign;
ChangeTwoLines(new, i, ln, i);
END;
FOR j := i+1 TO new.m-1 DO
skalar := new.elements[j, i]/pivot;
MulLineAddLine(new, i, j, i, skalar);
END;
END;
IF ~sign THEN
FOR i := 0 TO new.n-1 DO
new.elements[new.m-1, i] := -1.0 * new.elements[new.m-1, i];
END;
END;
RETURN new;
END Gauss;
PROCEDURE (a: MATRIX) Det * (): LONGREAL;
(*
* require
* a.DimM() = a.DimN()
*)
VAR
new: MATRIX;
i: LONGINT;
result: LONGREAL;
BEGIN
new := a.Gauss()(MATRIX);
result := new.elements[0, 0];
FOR i := 1 TO new.m-1 DO
result := result * new.elements[i, i];
END;
RETURN result;
END Det;
PROCEDURE (a: MATRIX) Norm * (): LONGREAL;
BEGIN
RETURN a.Det();
END Norm;
PROCEDURE (a: MATRIX) Compare * (b: BT.COMPAREABLE): LONGINT;
(*
* require
* a.DimM() = a.DimN()
* b(MATRIX).DimM() = b(MATRIX).DimN()
*)
VAR
i, j: LONGINT;
an, bn: LONGREAL;
BEGIN
WITH b: MATRIX DO
an := a.Norm(); bn := b.Norm();
IF an = bn THEN RETURN 0
ELSIF an > bn THEN RETURN 1
ELSE RETURN -1 END;
END;
END Compare;
END MATRIX.